#include "directsampling.h"
#include <cmath>

void direct_sampling_cont(float *retval, float * ti, int length, int width, int height, float *dv, int NR_variable, int NR_try, float tolerance, dv_bnd &bnd)
{
    for (int v_i=0;v_i<NR_variable;v_i++)
        retval[v_i]=0;
    int xsize=length+bnd.dvminx-bnd.dvmaxx;
    int ysize=width+bnd.dvminy-bnd.dvmaxy;
    int zsize=height+bnd.dvminz-bnd.dvmaxz;
    int dmax=255;
    float closest_diff=1024;
    for(int i=0;i<NR_try;i++){
        int xpos=rand()%xsize-bnd.dvminx;
        int ypos=rand()%ysize-bnd.dvminy;
        int zpos=rand()%zsize-bnd.dvminz;
        double diff=0;
        for(int dvi=0;dvi<dv[0];dvi++){
            for (int v_i=0;v_i<NR_variable;v_i++){
                diff+=fabs(ti[width*length*v_i+xpos+int(dv[1+dvi*(3+NR_variable)+0]) +length*(ypos+int(dv[1+dvi*(3+NR_variable)+1]))+length*width*(zpos+int(dv[1+dvi*(3+NR_variable)+2]))]-dv[1+(3+NR_variable)*dvi+3+v_i]);
            }
        }
        if(dv[0]>0)
            diff=diff/dmax/dv[0]/NR_variable;
        if (diff<=closest_diff) {
            for (int v_i=0;v_i<NR_variable;v_i++)
                retval[v_i]=ti[v_i*width*length+xpos+length*ypos+length*width*zpos];
            closest_diff=diff;
        }
        if (diff<tolerance) break;
    }
    return;
}

void Simulation_cont(Index <L2<float> > *R, float *ti, int NR_neighbor,int NR_try, float tolerance, vector<coord> &sample, vector<int> &ids, float *target, int N_sample, int tilength, int tiwidth, int tiheight, int NR_variable, int unitsize, int start_id, int level, int length, int width, int height)
{
    //The path is already random.
    float *dv=new float [1+(3+NR_variable)*NR_neighbor];
    //Data event [NR_cond, x,y,z,value, x,y,z,value, ...]
    for(int s_no=0;s_no<N_sample;s_no++){
        int fid=ids[s_no+start_id];//the real id in the sample
        vector <vector<int>> indices;
        getNearest_point(R,sample[fid].x,sample[fid].y,sample[fid].z,NR_neighbor, indices);
        size_t n_neigh=indices[0].size();
        if (n_neigh>NR_neighbor)
            dv[0]=NR_neighbor;
        else
            dv[0]=n_neigh;
        if ((s_no+start_id)<100) dv[0]=0;// At the beginning of unconditional simuation, set some cells to random values to avoid that the simulation result is too similar to the training image because of verbatim copy
        dv_bnd bnd;
        memset(&bnd,0,sizeof(int)*6);
        int dv_seq_no=0;
        for (int i=0;i<dv[0];i++) {
            int nfid=indices[0][i]-1;//the nfid shold be transformed to its real id, which is stored in ids
            nfid=ids[nfid];//get the real id
            if (nfid==fid) continue;
            //construct data event
            dv[1+dv_seq_no*(3+NR_variable)]=int((sample[nfid].x-sample[fid].x)/unitsize);
            dv[1+dv_seq_no*(3+NR_variable)+1]=int((sample[nfid].y-sample[fid].y)/unitsize);
            dv[1+dv_seq_no*(3+NR_variable)+2]=int((sample[nfid].z-sample[fid].z)/unitsize);
            for(int v_i=0;v_i<NR_variable;v_i++){
                dv[1+dv_seq_no*(3+NR_variable)+3+v_i]=target[v_i*length*width+sample[nfid].x+length*sample[nfid].y+length*width*sample[nfid].z];
            }
            if (dv[1+(3+NR_variable)*dv_seq_no]<bnd.dvminx) bnd.dvminx=dv[1+(3+NR_variable)*dv_seq_no];
            if (dv[1+(3+NR_variable)*dv_seq_no]>bnd.dvmaxx) bnd.dvmaxx=dv[1+(3+NR_variable)*dv_seq_no];
            if (dv[1+(3+NR_variable)*dv_seq_no+1]<bnd.dvminy) bnd.dvminy=dv[1+(3+NR_variable)*dv_seq_no+1];
            if (dv[1+(3+NR_variable)*dv_seq_no+1]>bnd.dvmaxy) bnd.dvmaxy=dv[1+(3+NR_variable)*dv_seq_no+1];
            if (dv[1+(3+NR_variable)*dv_seq_no+2]<bnd.dvminz) bnd.dvminz=dv[1+(3+NR_variable)*dv_seq_no+2];
            if (dv[1+(3+NR_variable)*dv_seq_no+2]>bnd.dvmaxz) bnd.dvmaxz=dv[1+(3+NR_variable)*dv_seq_no+2];
            dv_seq_no++;
            if (dv_seq_no>NR_neighbor) break;
        }
        float cur_val[NR_variable];
        direct_sampling_cont(cur_val,ti,tilength,tiwidth,tiheight,dv,NR_variable,NR_try,tolerance,bnd);
        for(int v_i=0;v_i<NR_variable;v_i++){
            target[v_i*length*width+sample[fid].x+length*sample[fid].y+length*width*sample[fid].z]=cur_val[v_i];
        }
        addPoint(R,sample[fid].x,sample[fid].y,sample[fid].z);
    }
    delete []dv;
}

void load_ti_3D(float ** ti, char *file, int *length, int *width, int *height)
{
    string type;
    int tmp;
    fstream fs;
    fs.open(file);
    fs>>*length;
    fs>>*width;
    fs>>*height;
    fs>>tmp;
    fs>>type;
    int size=(*length)*(*width)*(*height);
    *ti=new float [size];
    for(int i=0;i<size;i++){
        fs>>(*ti)[i];
    }
    fs.close();
}

void load_ti_tiff(float ** ti, char *file, int &length, int &width, int &NR_variable)
{
    //Read multi band tiff as the training image, limited to 2D case
    GDALDataset *poDS;
    GDALAllRegister();
    poDS = (GDALDataset *) GDALOpen( file, GA_ReadOnly );
    if( poDS == NULL ) {
        cout<<"Raster file open failed."<<endl;
        exit(1);
    }
    width=poDS->GetRasterXSize();
    length=poDS->GetRasterYSize();
    NR_variable=poDS->GetRasterCount();
    int size=width*length*NR_variable;
    *ti=new float [size];
    GDALRasterBand *poBand;
    for(int i=1;i<=NR_variable;i++){
        poBand=poDS->GetRasterBand(i);
        poBand->RasterIO( GF_Read, 0, 0, width, length, (*ti)+(i-1)*width*length, width, length, GDT_Float32, 0, 0 );
    }
    return;
}
